"""
    Radar

Observation operators for weather radar observations.

This module implements forward and adjoint operators for:
- Reflectivity (dBZ): Related to precipitation intensity
- Radial velocity (m/s): Doppler velocity along radar beam

# Physical Principles

## Reflectivity
Radar reflectivity depends on hydrometeor concentration and size:

    Z = Σ Dᵢ⁶

where D is particle diameter. The reflectivity factor in dBZ:

    dBZ = 10·log₁₀(Z/Z₀)

where Z₀ = 1 mm⁶/m³.

For rain:
    Z ≈ 200·R¹·⁶

where R is rainfall rate [mm/hr].

## Radial Velocity
Doppler velocity is the wind component along the radar beam:

    Vᵣ = u·sin(θ)·cos(φ) + v·cos(θ)·cos(φ) + w·sin(φ)

where:
- (u, v, w): wind components
- θ: azimuth angle from north
- φ: elevation angle from horizontal

# Forward and Adjoint Operators

Forward: (qrain, qsnow, u, v, w) → (Z, Vᵣ)
Adjoint: ∂J/∂(Z,Vᵣ) → ∂J/∂(qrain, qsnow, u, v, w)

The hydrometeor forward operator requires assumptions about
drop size distribution and particle characteristics.
"""
module Radar

using LinearAlgebra
using ..Interpolation

# Import interpolation functions
import ..Interpolation: trilinear_interpolate, trilinear_interpolate_adjoint!

export RadarOperator
export reflectivity_forward, reflectivity_adjoint
export radial_velocity_forward, radial_velocity_adjoint
export compute_reflectivity_from_rainwater, compute_reflectivity_from_snow
export dbz_to_linear, linear_to_dbz

# Physical constants
const rho_water = 1000.0    # Water density [kg/m³]
const rho_ice = 917.0       # Ice density [kg/m³]
const Z0 = 1.0              # Reference reflectivity [mm⁶/m³]
const pi = 3.141592653589793

"""
    RadarOperator{T<:AbstractFloat}

Weather radar observation operator configuration.

# Fields
- `observable::Symbol`: :reflectivity or :velocity
- `radar_lon::T`: Radar longitude [degrees]
- `radar_lat::T`: Radar latitude [degrees]
- `radar_alt::T`: Radar altitude [m]
- `beam_width::T`: Beam width [degrees]
- `wavelength::T`: Radar wavelength [m] (e.g., 0.05 for S-band)
"""
struct RadarOperator{T<:AbstractFloat}
    observable::Symbol
    radar_lon::T
    radar_lat::T
    radar_alt::T
    beam_width::T
    wavelength::T

    function RadarOperator{T}(;
                             observable::Symbol=:reflectivity,
                             radar_lon::T=T(0),
                             radar_lat::T=T(0),
                             radar_alt::T=T(0),
                             beam_width::T=T(1.0),
                             wavelength::T=T(0.1)) where T
        @assert observable in [:reflectivity, :velocity] "Invalid observable type"
        new{T}(observable, radar_lon, radar_lat, radar_alt, beam_width, wavelength)
    end
end

RadarOperator(; kwargs...) = RadarOperator{Float64}(; kwargs...)

# Reflectivity Operators

"""
    compute_reflectivity_from_rainwater(qr::T, rho_air::T) -> T

Compute radar reflectivity from rainwater mixing ratio.

Uses empirical relationship: Z = a·qrᵇ

# Arguments
- `qr`: Rainwater mixing ratio [kg/kg]
- `rho_air`: Air density [kg/m³]

# Returns
- Reflectivity factor Z [mm⁶/m³]
"""
function compute_reflectivity_from_rainwater(qr::T, rho_air::T) where T<:AbstractFloat
    if qr <= T(0)
        return T(0)
    end

    # Convert mixing ratio to rainwater content [g/m³]
    qr_content = qr * rho_air * T(1000)  # kg/m³ to g/m³

    # Empirical Z-qr relationship (simplified Marshall-Palmer)
    # Z = 3630 · qr^1.75 (qr in g/m³)
    Z = T(3630) * qr_content^T(1.75)

    return max(Z, T(0))
end

"""
    compute_reflectivity_from_snow(qs::T, rho_air::T) -> T

Compute radar reflectivity from snow mixing ratio.

Snow has different scattering properties than rain.
"""
function compute_reflectivity_from_snow(qs::T, rho_air::T) where T<:AbstractFloat
    if qs <= T(0)
        return T(0)
    end

    # Convert mixing ratio to snow content [g/m³]
    qs_content = qs * rho_air * T(1000)

    # Snow reflectivity (simplified)
    # Assumes lower reflectivity per unit mass than rain
    Z = T(2000) * qs_content^T(2.0)

    return max(Z, T(0))
end

"""
    linear_to_dbz(Z::T) -> T

Convert linear reflectivity to dBZ.

dBZ = 10·log₁₀(Z/Z₀)
"""
function linear_to_dbz(Z::T) where T<:AbstractFloat
    if Z <= T(0)
        return T(-30)  # Minimum detectable signal
    end

    return T(10) * log10(Z / Z0)
end

"""
    dbz_to_linear(dbz::T) -> T

Convert dBZ to linear reflectivity.

Z = Z₀·10^(dBZ/10)
"""
function dbz_to_linear(dbz::T) where T<:AbstractFloat
    return Z0 * T(10)^(dbz / T(10))
end

"""
    reflectivity_forward(qr_field::AbstractArray{T,3},
                        qs_field::AbstractArray{T,3},
                        rho_field::AbstractArray{T,3},
                        lon_grid::AbstractVector{T},
                        lat_grid::AbstractVector{T},
                        z_grid::AbstractVector{T},
                        lon_obs::T, lat_obs::T, z_obs::T) -> T

Forward operator for radar reflectivity observations.

# Arguments
- `qr_field`: Rainwater mixing ratio field [kg/kg]
- `qs_field`: Snow mixing ratio field [kg/kg]
- `rho_field`: Air density field [kg/m³]
- Grid coordinates and observation location

# Returns
- Reflectivity in dBZ
"""
function reflectivity_forward(qr_field::AbstractArray{T,3},
                              qs_field::AbstractArray{T,3},
                              rho_field::AbstractArray{T,3},
                              lon_grid::AbstractVector{T},
                              lat_grid::AbstractVector{T},
                              z_grid::AbstractVector{T},
                              lon_obs::T, lat_obs::T, z_obs::T) where T<:AbstractFloat

    # Interpolate fields to observation location
    qr_obs = trilinear_interpolate(qr_field, lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)
    qs_obs = trilinear_interpolate(qs_field, lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)
    rho_obs = trilinear_interpolate(rho_field, lon_grid, lat_grid, z_grid,
                                    lon_obs, lat_obs, z_obs)

    # Compute reflectivity contributions
    Z_rain = compute_reflectivity_from_rainwater(max(qr_obs, T(0)), rho_obs)
    Z_snow = compute_reflectivity_from_snow(max(qs_obs, T(0)), rho_obs)

    # Total reflectivity (sum of rain and snow contributions)
    Z_total = Z_rain + Z_snow

    # Convert to dBZ
    dbz = linear_to_dbz(Z_total)

    return dbz
end

"""
    reflectivity_adjoint!(qr_gradient::AbstractArray{T,3},
                         qs_gradient::AbstractArray{T,3},
                         rho_gradient::AbstractArray{T,3},
                         residual::T,
                         qr_field::AbstractArray{T,3},
                         qs_field::AbstractArray{T,3},
                         rho_field::AbstractArray{T,3},
                         lon_grid::AbstractVector{T},
                         lat_grid::AbstractVector{T},
                         z_grid::AbstractVector{T},
                         lon_obs::T, lat_obs::T, z_obs::T)

Adjoint operator for radar reflectivity observations.

Distributes reflectivity gradient back to hydrometeor fields.
"""
function reflectivity_adjoint!(qr_gradient::AbstractArray{T,3},
                              qs_gradient::AbstractArray{T,3},
                              rho_gradient::AbstractArray{T,3},
                              residual::T,
                              qr_field::AbstractArray{T,3},
                              qs_field::AbstractArray{T,3},
                              rho_field::AbstractArray{T,3},
                              lon_grid::AbstractVector{T},
                              lat_grid::AbstractVector{T},
                              z_grid::AbstractVector{T},
                              lon_obs::T, lat_obs::T, z_obs::T) where T<:AbstractFloat

    # Get interpolated values
    qr_obs = trilinear_interpolate(qr_field, lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)
    qs_obs = trilinear_interpolate(qs_field, lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)
    rho_obs = trilinear_interpolate(rho_field, lon_grid, lat_grid, z_grid,
                                    lon_obs, lat_obs, z_obs)

    # Compute reflectivities
    Z_rain = compute_reflectivity_from_rainwater(max(qr_obs, T(0)), rho_obs)
    Z_snow = compute_reflectivity_from_snow(max(qs_obs, T(0)), rho_obs)
    Z_total = Z_rain + Z_snow

    # Adjoint of dBZ conversion: ∂dBZ/∂Z = 10/(Z·ln(10))
    if Z_total > eps(T)
        ddbz_dZ = T(10) / (Z_total * log(T(10)))
    else
        ddbz_dZ = T(0)
    end

    residual_Z = residual * ddbz_dZ

    # Adjoint of reflectivity computations
    # ∂Z_rain/∂qr = 1.75 · 3630 · (qr·ρ·1000)^0.75 · ρ·1000
    if qr_obs > eps(T)
        qr_content = qr_obs * rho_obs * T(1000)
        dZ_dqr = T(1.75) * T(3630) * qr_content^T(0.75) * rho_obs * T(1000)
        qr_residual = residual_Z * dZ_dqr
    else
        qr_residual = T(0)
    end

    # ∂Z_snow/∂qs = 2.0 · 2000 · (qs·ρ·1000)^1.0 · ρ·1000
    if qs_obs > eps(T)
        qs_content = qs_obs * rho_obs * T(1000)
        dZ_dqs = T(2.0) * T(2000) * qs_content^T(1.0) * rho_obs * T(1000)
        qs_residual = residual_Z * dZ_dqs
    else
        qs_residual = T(0)
    end

    # Adjoint of trilinear interpolation
    trilinear_interpolate_adjoint!(qr_gradient, qr_residual,
                                   lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)
    trilinear_interpolate_adjoint!(qs_gradient, qs_residual,
                                   lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)

    # Density gradient (simplified - neglecting density dependencies)
    # Full implementation would include ∂Z/∂ρ terms

    return nothing
end

# Radial Velocity Operators

"""
    compute_radial_velocity(u::T, v::T, w::T,
                           azimuth::T, elevation::T) -> T

Compute radar radial velocity from wind components.

Vᵣ = u·sin(az)·cos(el) + v·cos(az)·cos(el) + w·sin(el)

# Arguments
- `u`, `v`, `w`: Wind components [m/s]
- `azimuth`: Azimuth angle from north [degrees]
- `elevation`: Elevation angle from horizontal [degrees]

# Returns
- Radial velocity [m/s] (positive = toward radar)
"""
function compute_radial_velocity(u::T, v::T, w::T,
                                azimuth::T, elevation::T) where T<:AbstractFloat

    # Convert angles to radians
    az_rad = deg2rad(azimuth)
    el_rad = deg2rad(elevation)

    # Compute radial velocity
    vr = u * sin(az_rad) * cos(el_rad) +
         v * cos(az_rad) * cos(el_rad) +
         w * sin(el_rad)

    return vr
end

"""
    radial_velocity_forward(u_field::AbstractArray{T,3},
                           v_field::AbstractArray{T,3},
                           w_field::AbstractArray{T,3},
                           lon_grid::AbstractVector{T},
                           lat_grid::AbstractVector{T},
                           z_grid::AbstractVector{T},
                           lon_obs::T, lat_obs::T, z_obs::T,
                           radar_lon::T, radar_lat::T, radar_alt::T) -> T

Forward operator for radar radial velocity observations.

# Arguments
- `u_field`, `v_field`, `w_field`: Wind component fields [m/s]
- Grid coordinates and observation location
- `radar_lon`, `radar_lat`, `radar_alt`: Radar location

# Returns
- Radial velocity [m/s]
"""
function radial_velocity_forward(u_field::AbstractArray{T,3},
                                v_field::AbstractArray{T,3},
                                w_field::AbstractArray{T,3},
                                lon_grid::AbstractVector{T},
                                lat_grid::AbstractVector{T},
                                z_grid::AbstractVector{T},
                                lon_obs::T, lat_obs::T, z_obs::T,
                                radar_lon::T, radar_lat::T, radar_alt::T) where T<:AbstractFloat

    # Interpolate wind components to observation location
    u_obs = trilinear_interpolate(u_field, lon_grid, lat_grid, z_grid,
                                  lon_obs, lat_obs, z_obs)
    v_obs = trilinear_interpolate(v_field, lon_grid, lat_grid, z_grid,
                                  lon_obs, lat_obs, z_obs)
    w_obs = trilinear_interpolate(w_field, lon_grid, lat_grid, z_grid,
                                  lon_obs, lat_obs, z_obs)

    # Compute beam geometry
    azimuth, elevation = compute_beam_geometry(radar_lon, radar_lat, radar_alt,
                                               lon_obs, lat_obs, z_obs)

    # Compute radial velocity
    vr = compute_radial_velocity(u_obs, v_obs, w_obs, azimuth, elevation)

    return vr
end

"""
    radial_velocity_adjoint!(u_gradient::AbstractArray{T,3},
                            v_gradient::AbstractArray{T,3},
                            w_gradient::AbstractArray{T,3},
                            residual::T,
                            lon_grid::AbstractVector{T},
                            lat_grid::AbstractVector{T},
                            z_grid::AbstractVector{T},
                            lon_obs::T, lat_obs::T, z_obs::T,
                            radar_lon::T, radar_lat::T, radar_alt::T)

Adjoint operator for radar radial velocity observations.

Distributes velocity gradient back to u, v, w wind fields.
"""
function radial_velocity_adjoint!(u_gradient::AbstractArray{T,3},
                                 v_gradient::AbstractArray{T,3},
                                 w_gradient::AbstractArray{T,3},
                                 residual::T,
                                 lon_grid::AbstractVector{T},
                                 lat_grid::AbstractVector{T},
                                 z_grid::AbstractVector{T},
                                 lon_obs::T, lat_obs::T, z_obs::T,
                                 radar_lon::T, radar_lat::T, radar_alt::T) where T<:AbstractFloat

    # Compute beam geometry
    azimuth, elevation = compute_beam_geometry(radar_lon, radar_lat, radar_alt,
                                               lon_obs, lat_obs, z_obs)

    # Convert angles to radians
    az_rad = deg2rad(azimuth)
    el_rad = deg2rad(elevation)

    # Adjoint of radial velocity computation
    # Vr = u·sin(az)·cos(el) + v·cos(az)·cos(el) + w·sin(el)
    # Therefore:
    # ∂Vr/∂u = sin(az)·cos(el)
    # ∂Vr/∂v = cos(az)·cos(el)
    # ∂Vr/∂w = sin(el)

    u_residual = residual * sin(az_rad) * cos(el_rad)
    v_residual = residual * cos(az_rad) * cos(el_rad)
    w_residual = residual * sin(el_rad)

    # Adjoint of trilinear interpolation
    trilinear_interpolate_adjoint!(u_gradient, u_residual,
                                   lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)
    trilinear_interpolate_adjoint!(v_gradient, v_residual,
                                   lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)
    trilinear_interpolate_adjoint!(w_gradient, w_residual,
                                   lon_grid, lat_grid, z_grid,
                                   lon_obs, lat_obs, z_obs)

    return nothing
end

# Beam Geometry Utilities

"""
    compute_beam_geometry(radar_lon::T, radar_lat::T, radar_alt::T,
                         obs_lon::T, obs_lat::T, obs_alt::T) -> Tuple{T,T}

Compute radar beam azimuth and elevation angles.

# Arguments
- Radar and observation locations

# Returns
- `(azimuth, elevation)`: Beam angles in degrees
"""
function compute_beam_geometry(radar_lon::T, radar_lat::T, radar_alt::T,
                              obs_lon::T, obs_lat::T, obs_alt::T) where T<:AbstractFloat

    # Simplified flat-earth geometry
    # For production, would use proper spherical geometry

    # Horizontal distance components (simplified)
    dx = (obs_lon - radar_lon) * T(111000) * cos(deg2rad(radar_lat))  # meters
    dy = (obs_lat - radar_lat) * T(111000)  # meters
    dz = obs_alt - radar_alt

    # Horizontal range
    range_horiz = sqrt(dx^2 + dy^2)

    # Azimuth from north (clockwise)
    azimuth = rad2deg(atan(dx, dy))  # atan2(x, y) for azimuth from north
    if azimuth < T(0)
        azimuth += T(360)
    end

    # Elevation angle
    elevation = rad2deg(atan(dz, range_horiz))

    return (azimuth, elevation)
end

"""
    compute_slant_range(radar_lon::T, radar_lat::T, radar_alt::T,
                       obs_lon::T, obs_lat::T, obs_alt::T) -> T

Compute slant range from radar to observation point.
"""
function compute_slant_range(radar_lon::T, radar_lat::T, radar_alt::T,
                            obs_lon::T, obs_lat::T, obs_alt::T) where T<:AbstractFloat

    # Simplified flat-earth
    dx = (obs_lon - radar_lon) * T(111000) * cos(deg2rad(radar_lat))
    dy = (obs_lat - radar_lat) * T(111000)
    dz = obs_alt - radar_alt

    return sqrt(dx^2 + dy^2 + dz^2)
end

end # module Radar
